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Abstract 



eigenvalues in [Eq, E\] is 



This paper presents numerical evidence that for quan- 
tum systems with chaotic classical dynamics, the num- 
ber of scattering resonances near an energy E scales like 

h 2 as h — > 0. Here, Kg denotes the subset of the 
classical energy surface {H = E} which stays bounded 
for all time under the flow generated by the Hamiltonian 
H and D (Ke) denotes its fractal dimension. Since the 
number of bound states in a quantum system with n de- 
grees of freedom scales like fr~ n , this suggests that the 



D(K E ) + 1 
2 

grees of freedom in scattering problems. 



quantity — ' - represents the effective number of de- 



1 Introduction 

Quantum mechanics identifies the energies of stationary 
states in an isolated physical system with the eigenvalues 
of its Hamiltonian operator. Because of this, eigenval- 
ues play a central role in the study of bound states, such 
as those describing the electronic structures of atoms and 
molecules.]] When the corresponding classical system al- 
lows escape to infinity, resonances replace eigenvalues as 
fundamental quantities: The presence of a resonance at 
A = E — vy, with E real and 7 > 0, gives rise to a dis- 
sipative metastable state with energy E and decay rate 7, 
as described in Such states are essential in scattering 
theory.0 

An important property of eigenvalues is that one can 
count them using only the classical Hamiltonian function 



H(x,p) 



V(x) and Planck's constant h: For 



fixed energies Eq < E\, the number N e i g (Eo, E\, K) of 
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1 For examples, see 

2 Systems which are not effectively isolated but interact only weakly 
with their environment can also exhibit resonant behavior. For example, 
electronic states of an "isolated" hydrogen atom are eigenfunctions of 
a self-adjoint operator, but coupling the electron to the radiation field 
turns those eigenstates into metastable states with finite lifetimes. This 
paper does not deal with dissipative systems and is only concerned with 
scattering. 



N eig {E ,E^h) 



vol({£o < H < Ei}) 



(1) 



where n denotes the number of degrees of freedom and 
vol(-) phase space volume. This result, known as the 
Weyl law, expresses the density of quantum states using 
the classical Hamiltonian function.]] No direct generaliza- 
tion to resonances is currently known. 

In this paper, numerical evidence for a Weyl-like power 
law is presented for resonances in a two-dimensional 
model with three symmetrically-placed gaussian poten- 
tials. A conjecture, based on the work of Sjostrand 
and Zworski [35|, states that the number of resonances 
A = E — ij with E < E < E\ and < 7 < Ti asymptot- 

+ l D ( K Eq) + 1 

ically lies between C\h 2 and C^h 2 
as h — > 0, where 



D (•) = dimension (see below), 

K E = K n {H = E} , 

K = {(x,p) : sup t ||$ t (x,p)|| < 00} , 

$4 = flow generated by H. 



(2) 



If D (Ke) depends continuously on E and \E\ — Eq\ is 
sufficiently small, then D (Ke x ) ~ D (Ke ) and the 
number of resonances in such a region is comparable to 

D(K E ) + 1 

h 2 for any E € [E , Ei\. 

The sets K and Ke are trapped sets and consist of 
initial conditions which generate trajectories that stay 
bounded forever. In systems where {H < E} is bounded 
for all E, the conjecture reduces to the Weyl asymptotic 

hr n . 

The notion of dimension requires some comment: The 
"triple gaussian" model considered here has very few 
trapped trajectories, and K and Ke (for any energy 
E) have vanishing Lebesgue measures. Thus, D(K) is 
strictly less than 2n — 4 and D (K E ) < 2n — 1 = 3. 
In fact, the sets K and Ke are fractal, as are trapped sets 
in many other chaotic scattering problems. Also, in this 
paper, the term "chaotic" always means hyperbolic; see 
Sjostrand [ 27 1 or Gaspard [ |l2| ] for definitions. 



3 For a beautiful exposition of early work on this and related themes, 
see [jl™. For recent work in the semiclassical context, see M. 
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This paper is organized as follows: First, the model sys- 
tem is defined. This is followed by mathematical back- 
ground information, as well as a heuristic argument for 
the conjecture. Then, numerical methods for computing 
resonances and fractal dimensions are developed, and nu- 
merical results are presented and compared with known 
theoretical predictions. 



Notation. 

function | | \p\ | 
tonian operator 



In this paper, H denotes the Hamiltonian 
' V(x) and H the corresponding Hamil- 

V{x), where A = Y^k- 



=1 as; 



the usual Laplacian and V acts by multiplication. 



2 Triple Gaussian Model 

The model system has n = 2 degrees of freedom; 
its phase space is i? 4 , whose points are denoted by 
(x,y,p x ,p y ). 
First, it is convenient to define 



G£ (*)=exp {-K^f 

Similarly, put 

GU,vo)( x 'V) = ( G *o®G° ya )(x,y) 

= G% o {x) ■ Gl Q {y) 

in two dimensions. 
Now, define H by 

H(x,y,p x ,p y ) = i (pI+P 2 v ) + V 3 (x,y) 

where the potential V m is given by 

V — V" r" 7 
c(k,m) = (iJcos(^) , J Rsin(^)). 



(3) 



(4) 



(5) 



(6) 



That is, it consists of rn gaussian "bumps" placed at the 
vertices of a regular m-gon centered at the origin, at a 
distance R > from the origin. This paper focuses 
on the case m = 3 because it is the simplest case that 
exhibits nontrivial dynamics in two dimensions. How- 
ever, the case m = 2 is also relevant because it is well- 
understood: See Mille r [|21| for early heuristic results and 
Gerard and Sjostrand [|13|| for a rigorous treatment. Thus, 
double gaussian scattering serves as a useful test case for 
the techniques described here. 

The quantized Hamiltonian H is similarly defined: 



H : 



(7) 




Figure 1 : Triple gaussian potential 

3 Background 

This section provides a general discussion of resonances 
and motivates the conjecture in the context of the triple 
gaussian model. However, the notation reflects the fact 
that most of the definitions and arguments here carry over 
to more general systems with n degrees of freedom. The 
reader should keep in mind that n = 2 for the triple gaus- 
sian model. 

There exists an extensive literature on resonances and 
semiclassical asymptotics in other settings. For example, 
see [|[ |o[ [TTJ |34| for detailed studies of the classical and 
quantum mechanics of hard disc scattering. 

3.1 Resonances 

Resonances can be defined mathematically as follows: Set 
R(z) = (H — zl)^ 1 for real z, where / is the identity op- 
erator. This one-parameter family of operators R(z) is the 
resolvent and is meromorphic with suitable modifications 
of its domain and range. The poles of its continuation into 
the complex plane are, by definition, the resonances of 

Less abstractly, resonances are generalized eigenval- 
ues of H. Thus, we should solve the time-independent 
Schrodinger equation 



Htp = Xip 



(8) 



to obtain the resonance A and its generalized eigenfunc- 
tion ip. In bound state computations, one approximates 
^ as a finite linear combination of basis functions and 
solves a finite-dimensional version of the equation above. 



See Figure 111. 



4 For more details and some references, see Jj^ ] . 



2 



To carry out similar calculations for resonances, it is nec- 
essary that ip lie in a function space which facilitates such 
approximations, for example L 2 . 

Let ip and A solve (||). Then e~^ At • ip solves the time- 
dependent Schrodinger equation 



in 9 4 = m. 

at 



(9) 



It follows that Im (A) must be negative because 
metastable states decay in time. Now suppose, for sim- 
plicity, that n = In Then solutions of (M) with energy 



E behave like e s 



/2Ex 



for large x > 0. Substituting 



A = E — ij for E yields e * 
ponentially because Im (V-^ 1 — 



2A:c 



which grows ex- 
< 0. Thus, finite 



rank approximations of H cannot capture such general- 
ized eigenfunctions. However, if we make the formal 
substitution x h—s- xe la , then the wave function becomes 



exp 



2A- 



Choosing a > |tan 1 (j/E) 
forces ip to decay exponentially. 

This procedure, called complex scaling, transforms the 
Hamiltonian operator H into the scaled operator H a . 
It also maps metastable states ip with decay rate 7 < 
E tan (2a) to genuine L 2 eigenfunctions ip a of H a . The 
corresponding resonance A becomes a genuine eigen- 
value: H a ip a = \ipa- Furthermore, resonances of H 
will be invariant under small perturbations in a, whereas 
other eigenvalues of H a will not. The condition a > 
\ tan -1 ('y/E) implies that, for small 7 and fixed E, the 
method will capture a resonance A = E — 17 if and only 
if 7 < 2Ea + O (a 2 ). We can perform complex scaling 



in polar where 



in higher dimensions by substituting r 1— > r 
coordinates. 

In algorithmic terms, this means we can compute eigen- 
values of H a for a few different values of a and look for 
invariant values, as demonstrated in Figure ||. In addition 
to its accuracy and flexibility, this is one of the advan- 
tages of complex scaling: The invariance of resonances 
under perturbations in a provides an easy way to check 
the accuracy of calculations, mitigating some of the un- 
certainties inherent in computational work.[] Note that the 
scaled operator H a is no longer self-adjoint, which results 
in non-hermitian finite-rank approximations and complex 
eigenvalues. 

This method, first introduced for theoretical purposes 
by Aguilar and Combes [[[J and Balslev and Combes [j^], 
was further developed by B. Simon in [p6[|. It has since 
become one of the main tools for computing resonances 
in physical chemistry |p2|, [HL p2[ g4j ]. For recent mathe- 
matical progress, see |]1 & [27], p8| and references therein. 

5 The analysis in higher dimensions requires some care, but the es- 
sential result is the same. 

6 For a different approach to computing resonances, see |R3[ and the 
references there. 
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Figure 2: Illustration of complex scaling: The three lines 
indicate the location of the rotated continuous spectrum 
for different values of a, while the box at the top of the fig- 
ure is the region in which resonances are counted. Eigen- 
values which belong to different values of a are marked 
with different styles of points. As explained later, only 
eigenvalues near the region of interest are computed. This 
results in a seemingly empty plot. 



For reference, the scaled triple-gaussian operator H a is 



H a = -e- 2ia -—A + V 3 , a , (10) 



" m.OL 

c a (k,m) = e~ ia c(k,m) 



(ID 



Note that these expressions only make sense because 
G% Q (x) is analytic in x, xq, and a. 

3.2 Fractal Dimension 

Recall that the Minkowski dimension of a given set U C 
i? m is 

D = inf jd : limsup (e d ~ rn ■ vol (U £ )) < 00 j , (12) 

where U e — {y £ R m : dist(y, U) < e}. A simple calcu- 
lation yields 

D(U) = lim A 7 ") (13) 



*o log (1/e) 



if the limit exists. 
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Texts on the theory of dimensions typically begin 
with the Hausdorff dimension because it has many de- 
sirable properties. In contrast, the Minkowski dimension 
can be somewhat awkward: For example, a countable 
union of zero-dimensional sets (points) can have positive 
Minkowski dimension. But, the Minkowski dimension is 
sometimes easier to manipulate and almost always easier 
to compute. It also arises in the heuristic argument given 
below. 

For a detailed treatment of different definitions of di- 
mension and their applications in the study of dynamical 
systems, see M, Oj. 



3.3 Generalizing the Weyl Law 

The formula 

vol{{E a < H < Ei}) 



{2irK) n 



makes no sense in scattering problems because the vol- 
ume on the right hand side is infinite for most choices of 
Eq and Ei, and this seems to mean that there is no gener- 
alization of the Weyl law in the setting of scattering the- 
ory. However, the following heuristic argument suggests 
otherwise: 

As mentioned before, a metastable state corresponding 
to a resonance A = E — i-f has a time-dependent fac- 
tor of the form e~% xt — e~% Et ■ e~n*. A wave packet 
whose dynamics is dominated by A (and other resonances 
near it) would therefore exhibit temporal oscillations of 
frequency 0(E/T%) and lifetime 0(h/^). Heuristically, 
then, the number of times the particle "bounces" in the 
"trapping region']] before escaping should be comparable 

n 7 7 

In the semiclassical limit, the dynamics of the wave 
packet should be well-approximated by a classical trajec- 
tory. Let T(x,y,p x ,p y ) denote the time for the parti- 
cle to escape the system starting at position (x, y) with 
momentum (p x ,p y ). The diameter of the trapping re- 
gion is 0(R), and typical velocities in the energy surface 
{H = E} are 0(^/E) (mass set to unity), so the num- 
ber of times a classical particle bounces before escaping 
should be 0(T\f~E j 'R). This suggests that, in the limit 
7i — > 0, T^/E/R r-j E/"/ and consequently 



T 



rVe 

1 



(15) 



Fix 70 > 0, and consider 

Nres = #{E-ij:E <E KE^jKja} (16) 



7 For our triple gaussian system, that would be the triangular region 
bounded by the gaussian bumps. 



for fixed energies Eq and E\. Equation ( |15[ ) implies that 
T > R^/Eq/jq, so by analogy with the Weyl law, 



vol({£ <H<E U T> ^}) 



(17) 



follows as an approximation for the number of quantum 
states with the specified energies and decay rates. 

Now, the function 1/T is nonnegative for all 
(x,y,p x ,p y ) and vanishes on K[ Eo<El ] = K n 
{^o < H < Ei}. Assuming that 1/T is sufficiently reg- 
ular,^] this suggests 

1/T(x,y,p x ,p y ) ~ d K[Eo Ei] {x,y,p x ,p y ) 2 , (18) 

where dx [E B] denotes distance to K[ EgEl y It follows 
that N res should scale like 



(14) vol({^o<g<^i,^ [Eo ^ <7p*}) 



For small 70, this becomes 



2n-D K, 



C(R,E ,Ei)-h- n - 7o 



'll) 



(19) 



(20) 



for some constant C, by (jj~2|). Choosing 70 = h and 
assuming that D (K E ) decreases monotonically with in- 
creasing E (as is the case in Figure ^2|), we obtain 

o(k Ei )+i d(k Bo )+i 

Cih 5 — < N res < c h 5 — . (21) 

If \Ei — E \ is sufficiently small, then D (Ky Eo « 
D (K E ) + 1 for E G [E Q , Ei], and 



(22) 



In [p7[], Sjostrand proved the following rigorous upper 
bound: For 70 > satisfying Cfi < 70 < 1/C, 

N res = O lc s h- n lo 2 (23) 

holds for all S > 0. When the trapped set is of pure 
dimension, that is when the infimum in Equation ( |T2"| ) is 
achieved, one can take 5 = 0. Setting 70 = h gives an 
upper bound of the form (j22|). 

8 In fact, this is numerically self-consistent: Assume that 1/T van- 
ishes to order u (with v not necessarily equal to 2) on K, and as- 
sume the conjecture. Then the number of resonances would scale like 
fo(2n-D(K ))/i^ f rom which one can solve for v. With the numerical 
data we have, this indeed turns out to be 2 (but with significant fluctua- 
tions). 

Also, if 1/T does not vanish quadratically everywhere on K, varia- 
tions in its regularity may affect the correspondence between classical 
trapping and the distribution of resonances. 
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In his proof, Sjostrand used the semiclassical argument 
above with escape functions and the Weyl inequality for 
singular values. Zworski continued this work in []3~5j], 
where he proved a similar result for scattering on convex 
co-compact hyperbolic surfaces with no cusps. His work 
was motivated by the availability of a large class of ex- 
amples with hyperbolic flows, easily computable dimen- 
sions, and the hope that the Selberg trace formula could 
help obtain lower bounds. But, these hopes remain unful- 
filled so far Jl4|], and that partly motivates this work. 

4 Computing Resonances 

Complex scaling reduces the problem of calculating res- 
onances to one of computing eigenvalues. What remains 
is to approximate the operator H a by a rank N operator 
Hn,c and to develop appropriate numerical methods. For 
comparison, see [[22| [3l| |24j] for applications of com- 
plex scaling to problems in physical chemistry. 

4.1 Choice of Scaling Angle. 

One important consideration in resonance computation is 
the choice of the scaling angle a: Since we are inter- 
ested in counting resonances in a box [Eo, Ei] — i [0, fi], 

it is necessary to choose a > tan -1 (^j^j so that the 

continuous spectrum of H a is shifted out of the box 
[£ ,£i]-«[0,?i] (see Figure |). 
In fact, the resonance calculation uses 



= tan" 1 



(24) 



This choice of a helps avoid the pseudospectrum [30 36| ]: 
Let A be an N x N matrix, and let R(z) be the resolvent 
(A — zl)~ x . It is well known that when A is normal, 
that is when A commutes with its adjoint A*, the spectral 
theorem applies and the inequality 



\\R(*)\\ 



< 



11(4- zl)- 1 
dist(z,a(A)Y 



(25) 



holds (<j(A) denotes the spectrum of A). When A is not 
normal, no such inequality holds and | |i?(.z)| | can become 
very large for z far from <r(A). This leads one to define 
e-pseudo spectrum: 



K(A) = {z:\\R(z)\\>l/e}. 



(26) 



Using the fact that A is a matrix, one can show that A e (A) 
is equal to the set 

{z : 3A' such that z e a(A + A'),\\A'\\ < e} . (27) 



That is, the e-pseudospectrum of A consists of those com- 
plex numbers z which are eigenvalues of an e-perturbation 
of A. 

The idea of pseudospectrum can be extended to general 
linear operators. In p0|], it is emphasized that for non- 
normal operators, the pseudospectrum can create "false 
eigenvalues" which make the accurate numerical compu- 
tation of eigenvalues difficult. In [|36j], this phenomenon 
is explained using semiclassical asymptotics. Roughly 
speaking, the pseudospectrum of the scaled operator H a 
is given by the closure of 



{z : z = H a (x,y,p x ,p y )} 



(28) 



of its symbol H a , which is the scaled Hamiltonian func- 
tion 



H a (x,y,p x ,p y ) 



-2ia Px 



P 



V -L. V 

r V3,q 



(29) 



in this case. Choosing a to be comparable to fi ensures 
that the imaginary part of H a is also comparable to fi, 
which keeps the pseudospectrum away from the count- 
ing box [Eo, Ei] — i [0, h]; a larger a would contribute a 
larger a 2 term to the imaginary part of H a and enlarge 
the pseudospectrum. As one can see in Figures [3j] - [34| 
the invariance of resonances under perturbations in a also 
helps filter out pseudospectral effects. 

This consideration also points out the necessity of the 
choice a = tan -1 (tJ^) 1 To avoid pseudospectral ef- 
fects, a must be O (fi). On the other hand, if a — o (fi), 
then finite rank approximations may fail to capture reso- 
nances in the region of interest. 

4.2 Eigenvalue Computation 

Suppose that we have constructed H^, a . In the case of 
eigenvalues, the Weyl law states that N e i g = O (fi ) as 
fi ^ 0, since our system has n = 2 degrees of freedom. 
Thus, in order to capture a sufficient number of eigenval- 
ues, the rank N of the matrix approximation must scale 
like fC 2 . In the absence of more detailed information on 
the density of resonances, the resonance computation re- 
quires a similar assumption to ensure sufficient numerical 
resolution. 

Thus, for moderately small fi, the matrix has iV 2 ~ Ti~ A 
entries, which rapidly becomes prohibitive on most com- 
puters available today. Furthermore, even if one does not 
store the entire matrix, numerical packages like LAPACK 
[^] require 0(N 2 ) auxiliary storage, again making practi- 
cal calculations impossible. 

Instead of solving the eigenvalue problem H n a v = Xv 
directly, one solves the equivalent eigenvalue problem 



= \'v. 



(30) 
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Efficient implementations of the Arnoldi algorithm 
[ |T^ 1 can solve for the largest few eigenvalues A' of 

\Htf, a — -V)) ■ But A = Aq + 1/A', so this method 



allows one to compute a subset of the spectrum of H^ a 
near a given Ao. 

Such algorithms require a method for applying the ma- 
trix (-ffjv.a — to a given vector v at each iteration 
step. In the resonance computation, this is done by solv- 
ing (HN,a — Ao)u> = v for w by applying conjugate gra- 
dient to the normal equations (see [@])|] The resonance 
program, therefore, consists of two nested iterative meth- 
ods: An outer Arnoldi loop and an inner iterative linear 
solver for (H^ a — \o)w = v. This computation uses 
ARPACKp], which provides a flexible and efficient im- 
plementation of the Arnoldi method. 

To compute resonances near a given energy E, the pro- 
gram uses Ao = E + ia, a > 0, instead of Ao = E: This 
helps control the condition number of Hpj !a —Ao and gives 
better error estimates and convergence criteria]^] 

4.3 Matrix Representations 
4.3.1 Choice of Basis 

While one can discretize the differential operator H a 
via finite differences, in practice it is better to represent 
the operator using a basis for a subspace of the Hilbert 
space L 2 : This should better represent the properties of 
wave functions near infinity and obtain smaller (but more 
dense) matrices. 

Common basis choices in the chemical literature in- 
clude so-called "phase space gaussian" [^] and "dis- 
tributed gaussian" bases fll5|]. These bases are not orthog- 
onal with respect to the usual L 2 inner product, so one 
must explicitly orthonormalize the basis before comput- 
ing the matrix representation of H a . In addition to the 
computational cost, this also requires storing the entire 
matrix and severely limits the size of the problem one can 



'That is, instead of solving Aw = v, one solves A* Aw = A* v. 
This is necessary because -ffjv, Q is non-hermitian, and conjugate gradi- 
ent only works for positive definite matrices. This is not the best numer- 
ical method for non-hermitian problems, but it is easy to implement and 
suffices in this case. 

10 See for details on the package, as well as an overview of Krylov 
subspace methods. 

"Most of the error in solving the matrix equation (-ffjv, Q — Ao)ui = v 

concentrates on eigenspaces of (-Hjv a — -V)) _1 with large eigenvalues. 
These are precisely the desired eigenvalues, so in principle one can toler- 
ate inaccurate solutions. However, the calculation requires convergence 
criteria and error estimates for the linear solver, and using a > 0, say 
a = 1, turns out to ensure a relative error of about 10 — 6 after about 
17-20 iterations of the conjugate gradient solver. Since we only wanted 
to count eigenvalues, a more accurate (and expensive) computation of 
resonances was not necessary. 




Figure 3: A sine function with to = 0, Ax = 1. 



solve. Instead, this computation uses a discrete-variable 
representation (DVR) basis j20|]: 

Consider, for the moment, the one dimensional problem 
of finding a basis for a "good" subspace of L 2 (R). Fix a 
constant Ax > 0, and for each integer to, define 



sin (-r-(x — to Ax)) 
'Ax ^ . > " . (31) 



7r(x — toAx) 



(This is known as a "sine" function in engineering litera- 
ture [S]. See Figure |[) The Fourier transform of <fi m ,Ax 
is 



»m.Ai W) 



/ Ax, \w\ < ir/Ax 
\lu\ > it/ Ax 



(32) 



One can easily verify that {<fi m .Ax} forms an orthonor- 
mal basis for the closed subspace of L 2 functions whose 
Fourier transforms are supported in [— ir/Ax, ir/Ax]. 

To find a basis for corresponding space of band-limited 
functions in L 2 (R 2 ), simply form the tensor products 



(33) 



The basis has a natural one-to-one correspondence with 
points (toAi + Xq, 71A2 + Yq) on a regular lattice of grid 
points in a box [Xq,Xi] x [Yqj Y"i] covering the spatial 
region of interest. (See Figure ^.) Using this basis, it is 
easy to compute matrix elements for H a . 

4.3.2 Tensor Product Structure 

An additional improvement comes from the separability 
of the Hamiltonian: Each term in the scaled Hamiltonian 
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Figure 4: Illustration of resonance program parameters in 
configuration space: The lower-left corner of the mesh is 
(Xq, Yq), while the upper right corner is (X%, Y\). The 
mesh contains N x x N y grid points, and a basis function 
<t>mn is placed at each grid point. Stars mark the centers of 
the potentials, the circles have radius 2a (with a = 1/3), 
and R is set to 1.4. Parameters for the classical computa- 
tion are depicted in Figure 



H a splits into a tensor product: 



G 



dx 2 

dy 2 

(xo,yo) 



d2 r 
r d2 

dy 2 



Ct 



G° 



(34) 

(35) 
(36) 



where I x and I y denote identity operators on copies of 
L 2 (R). Since the basis {</>,„„} consists of tensor products 
of one dimensional bases, iJjv.a is also a short sum of 
tensor products. Thus, if we let N x denote the number of 
grid points in the x direction and let N y denote the number 
of grid points in the y direction, then N = 
Hjs[. a is a sum of five matrices of the form A 



N x ■ N y and 
; Cg> A v , where 



A x is N x x N x and A y is N y x N y . 

Such tensor products of matrices can be applied to ar- 
bitrary vectors efficiently using^ the outer product repre- 
sentation.^] Since the rank of Hn , a is TV = N x ■ N y and 
N x sa N y in these computations, we can store the tensor 
factors of the matrix i/jv,a using O(N) storage instead 
of 0{N 2 ), and apply H^.a to a vector in O (iV 3 / 2 ) time 
instead of O (N 3 ). The resulting matrix is not sparse, as 



12 The tensor product of two column vectors v and w can be repre- 
sented as v ■ w T . We then have (A 45 B) ■ (v J4 to) (Av) ■ (Bw) T , 



one can see from the matrix elements for the Laplacian 
below. 

Note that this basis fails to take advantage of the dis- 
crete rotational symmetry of the triple gaussian Hamil- 
tonian. Nevertheless, the tensor decomposition provides 
sufficient compression of information to facilitate efficient 
computation. 

4.3.3 Matrix Elements 

It is straightforward to calculate matrix elements for the 
Laplacian on R 1 : 



K„ 



Ax' 2 ' 



Ax 2 -(m — n) 2 ' 



m = n 
to n 



(37) 



There is no closed form expression for the matrix ele- 
ments of the potential, but it is easy to perform numerical 
quadrature with these functions. For example, to compute 



V„. 



G(x)4>m{x)4> n {x) dx 



(38) 



for G(x) = e 3^, one computes 

N 



V,r: 



J2 G(kS) ■ S ■ cf> m {k5) ■ <t> n (kS), (39) 



k=-N 



where the stepsize 5 should satisfy <5 < Ax/2. It is easy 
to show that the error is bounded by the sum of 



2 exp 



2S 2 



(40) 



which extends by linearity to (A £g) B) ■ u = A ■ 



B T . 



which controls the aliasing error, and 

2 ^M 2 / (iV-l)V\ 

. exp ^ . (41) 

Ax l \ 2\a\ 2 J 

which controls the truncation error. 

4.3.4 Other Program Parameters 

The grid spacing Ax implies a limit on the maximum pos- 
sible momentum in a wave packet formed by this basis. In 
order to obtain a finite-rank operator, it is also necessary 
to limit the number of basis functions. 

The resonance computation used the following param- 
eters: 

1. Xq, X±, Yq, and Y\ are chosen to cover the region of 
the configuration space for which Vj,(x, y) > 10~ 4 . 
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2. Let L x = X\ — Xq and L y = Y\ — Yq denote the 
dimension of the computational domain. The reso- 
nance calculation uses TV = N x ■ N y basis functions, 

L y VSE 



with AT T = 1.6 



2irh 



an&Ny = 1.6 



2nh 



3. This gives 

Ax = L x /N x , 
Ay = L y /N v , 

which limits the maximum momentum in a wave 



(42) 



packet to \p x \ 
irh/Ay = 1.6 



< irh/Ax = 1.6a 
[2E. 



zE and \p y \ < 



5 Trapped Set Structure 

5.1 Poincare Section 

Because the phase space for the triple gaussian model is 
R A and its flow is chaotic, a direct computation of the 
trapped set dimension is difficult. Instead, we try to com- 
pute its intersection with a Poincare section: 

Let E be a fixed energy, and recall that R is the distance 
from each gaussian bump to the origin. Choose R < 
R so that the circles Ck of radius R centered at each 
potential, for k = 0, 1,2, do not intersect. The angular 
momentum pg with respect to the /cth potential center is 
defined by pg = Ax ■ p y — Ay • p x , where Aa; = x — 
Rcos(9 k ), Ay = y- R$in(6 k ), and k = 

Let P be the submanifold P U Pi U P 2 of R 4 (see 
Figure |]), where the coordinates (0,pe) in the subman- 
ifold P k are related to ambient phase space coordinates 
(x,y,p x ,p y ) by 



x 

y 

Px 

Py 



R cos I 
R sin ( 

p r COS 

p r sin 



Ok) + Ro cos( 
9 k ) + R sin( 



k ), 

Ok) , 



9 k )-^-sm(9 + 6 k ) 



Ro 

> + e k ) + j£;Cos(8 + 9 k ) 



(43) 



and the radial momentum p r is 



Pr 



E - Vs{x ' y) ~wl' 



(44) 




-1 1 



Figure 5: A typical trajectory: Stars mark the potential 
centers. In this case, R — 1.4 and E = 0.5. The circles 
drawn in the figure have radius 1, and the disjoint union 
of their cotangent bundles form the Poincare section. Tra- 
jectories start on the circle centered at bump (the bumps 
are, counterclockwise, 0, 1, and 2) with some given angle 
9 and angular momentum pg. This trajectory generates the 
finite sequence (0, 1,2,0, 2, oo). (Symbolic sequences are 
discussed later in the paper.) An illustration of resonance 
computation is depicted in Figure The dashed line is the 
time-reversed trajectory with the same initial conditions, 
generating the sequence (oo, 2, 0, 2, 0). 



Note that this implicitly embeds P into the energy surface 
{H = E}, and the radial momentum p r is always posi- 
tive: The vector (p x ,p y ) points away from the center of 
C k - 

The trapped set is naturally partitioned into two sub- 
sets: The first consists of trajectories which visit all three 
bumps, the second of trajectories which bounce between 
two bumps. The second set forms a one-dimensional sub- 
space of Ke, so the finite stability of the Minkowski di- 
mension]^ implies that the second set does not contribute 

13 That is, D(A UB) = max (D(A), D(B)). For details, see [EJ|. 
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to the dimension of the trapped set. More importantly, 
most trajectories which visit all three bumps will also cut 
through P. 

One can thus reduce the dimension of the problem by 
restricting the flow to Ke H P, as follows: Take any 
point (9,pe) in Pk, and form the corresponding point 
(x,y,p x ,p y ) in i? 4 via Equation (j43|). Follow along the 
trajectory $ t (x,y,p x ,p y ). If the trajectory does not es- 
cape, eventually it must encounter one of the other circles, 
say Ck< ■ Generically, trajectories cross Cy twice at each 
encounter, and we denote the coordinates (0',p' e ) (in Py) 
of the outgoing intersection by 



,p e ,k) = (9',p' e ,k'). 



(45) 



If a trajectory escapes from the trapping^ region, we can 
symbolically assign oo to $. The map $ then generates 
stroboscopic recordings of the flow $ f on the submani- 
fold P, and the corresponding discrete dynamical system 
has trapped set Ke H P. So, instead of computing $ t on 
R 4 , one only needs to compute $ on P. By symmetry, it 
will suffice to compute the dimension of Ke = Ke D Pq. 
Pushing Ke forward along the flow $f adds one dimen- 
sion, so D (Ke) = D I Ke ) + 1. Being a subset of the 



two-dimensional space Pq, Ke is easier to work with. 

Readers interested in a more detailed discussion of 
Poincare sections and their use in dynamics are referred 
to |^9|]. For an application to the similar but simpler set- 
ting of hard disc scattering, see fyh. Also, Knauf has 
applied some of these ideas in a theoretical investigation 
of classical scattering by Coulombic potentials fll7|], 

5.2 Self-Similarity 

Much is known about the self-similar structure of the 
trapped set for hard disc scattering O]; less is known 
about "soft scatterers" like the triple gaussian system. 
However, computational results and analogy with hard 
disc scattering give strong support to the idea that K (and 
hence Ke) is self-similar.|^ Consider Figures ^ - [l2| : 
They show clearly that Ke is self-similar. (In these im- 
ages, E — 0.5 and Rq = 1-0.) However, it is also clear 
that, unlike objects such as the Cantor set or the Sierpinski 
gasket, Ke is not exactly self-similar. 




-0.5 



Figure 6: Points in Pq which do not go to 00 after one 
iteration of <£>. The horizontal axis is 9 and the vertical 
axis is pq. 




-0.5 



5.3 Symbolic Dynamics 

The computation of D (^Ke^J uses symbolic sequences, 
which requires a brief explanation: For any point (9,pg), 
let Si denote the third component of <& % (9,pe) (see (p3|)), 
for any integer i. Thus, Si is the index k of the circle 



4 More precisely, self-affine. 



Figure 7: Points in Pq which do not go to 00 after one 
iteration of The horizontal axis is 9 and the vertical 
axis is pg. 
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3.55 3.6 3.65 3.7 



Figure 8: The intersection of the sets in Figures || and Figure 10: The cut-out part of Figure |9| magnified. Recall 
0. These points correspond to symmetric sequences of that these correspond to symmetric sequences of length 5; 
length 3. compare with Figure [sj. 
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3.2 3.4 3.6 3.8 4 



Figure 9: The lower-right "island" in Figure |L magnified. Figure 1 1 : The upper-right island in Figure n. The white 
The white cut-out in the middle is the subset correspond- cut-out in the middle is, again, the subset corresponding 
ing to symmetric sequences of length 5. to symmetric sequences of length 5. 
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0.15 
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3.55 
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Figure 12: The cut-out part of Figure |TT[ magnified. 
Recall that these correspond to symmetric sequences of 
length 5. Compare with Figures Hand HQ. 



Cfe that the trajectory intersects at the ith iteration of $ 
(or the |i|th iteration of Such symbolic sequences 

s = (...Si, s i+1 ...) satisfy Si £ {0, 1, 2, oo} and Si ^ s i+ i 
for all i, and with oo occuring only at the ends. Let us call 
sequences satisfying these conditions valid. 

For example, the trajectory in Figure ^| generates the 
valid sequence (0, 1, 2, 0, 2, oo), where the dot over in- 
dicates that the initial point (6,pg) of the trajectory be- 
longs to Pq. Thus, we can label collections of trajectories 
using valid sequences, and label points in P with "dot- 
ted" sequences. Clearly, trapped trajectories generate bi- 
infinite sequences.^ 

The islands in Figures |] - [FT] correspond to sym- 
metric sequences centered at 0, of the form s = 
(...s-k, s_i,6, si, s 2 , s fc ...): By keeping track of 
the symbolic sequences generated by each trajectory, one 
can easily label and isolate each island. This is a useful 
property from the computational point of view. 

5.4 Dimension Estimates 

To compute the Minkowski dimension using Equation 
(pj|), jwe need to determine when a given point is within 
e of K E . This is generally impossible: The best one can 
do is to generate longer and longer trajectories which stay 
trapped for increasing (but finite) amounts of time. 

Instead, one can estimate a closely related quantity, the 
information dimension, in the following way: Let K L 



15 In hard disc scattering, the converse holds for sufficiently large R: 
To each bi-infinite valid sequence there exists a trapped trajectory gener- 
ating that sequence. This may not hold in the triple gaussian model, and 
in any case it is not necessary for the computation. 



denote the set of all points in Pq corresponding to sym- 
metric sequences of length 2k + 1 centered at 0. That is, 

~ (k) 

K E ' consists of all points in Pq which generate trajec- 
tories (both forwards and backwards in time) that bounce 

~ (k) 

at least k times before escaping. The sets K E decrease 
monotonically to K E : D K [ E +1) and rff =0 K% ) = 
K E . 

One can then estimate the information dimension using 
the following algorithm: 

1. Initialization: Cover Pq with a mesh Lq with Nq x 
No grid points and mesh size eo- 

2. Recursion: Begin with K E , which consists of four 
islands corresponding to symmetric sequences of 
length 2-1 + 1 = 3 (see Figure ||). Magnify each 
of these islands and compute the sub-islands corre- 
sponding to symmetric sequences of length 5 (see 
Figures ^ and [ll]). Repeat this procedure to recur- 
sively compute the islands of K E from those of 

fco, where ko is sufficiently 
large that each island of K E °^ has diameter smaller 
than the mesh size e of Lq. 



' (k) 

K E . Continue until k 



Estimation: Using the islands of K E 
probability 



(fco) 



vol [K { E a) C\B., 



vol 



estimate the 



(46) 



for the (ij')th cell of Lq. 
dimension via 



We can then compute the 



D 



(*■) 



(47) 



which reduces to 
form because eo ^ 



log(iVo) 
when the distribution is uni 



1/JVo 



Under suitable conditions (as is assumed to be the case 
here), the information dimension agrees with both the 
Hausdorff and the Minkowski dimensions.^ 

The algorithm begins with the lattice Lq with which 
one wishes to compute the dimension. It then recursively 

~ (k) 

computes K E for for increasing values of k, until it 
closely approximates Kb relative to the mesh size of Lq. 
It is easy to keep track of points belonging to each is- 
land in this computation, since each island corresponds 
uniquely to a finite symmetric sequence. Note that while 
the large mesh Lq remains fixed throughout the computa- 
tion, the recursive steps require smaller N\ x N\ meshes 

~ (k) 

around each island of K E up to the value of k specified 
by the algorithm. See Figure |l3[ 

16 See ||23|] for a discussion of the relationship between these dimen- 
sions, as well as their use in multifractal theory. 
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hbar=0.025 



Figure 13: This figure illustrates the recursive step in the 
dimension estimation algorithm: The dashed lines repre- 
sent Lq, while the solid lines represent a smaller mesh 
centered on one of the islands. The No x No mesh Lo 
remains fixed throughout the computation, but the smaller 
Ni x Ni mesh is constructed for each island of Kg up 
to the value of k specified by the algorithm. 
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Figure 14: These are the eigenvalues of H^^a, for 
E = 0.5, h = 0.025, R = 1.4, and a € 
{0.0624, 0.0799, 0.0973}. This calculation used an 102 x 
108 grid, and 90 out of N = 11016 eigenvalues were 
computed. 



6 Numerical Results 



6.1 Resonance Counting 



As an illustration of complex scaling, Figures [Ij - [18] 
contain resonances for R — 1.4 and h G [0.017,0.025]. 
Eigenvalues of H^.a for different values of a are marked 
by different styles of points, and the box has depth h and 
width 0.2, with E = 0.4 and E Y = 0.6. These plots may 
seem somewhat empty because only those eigenvalues of 
-Hjv a m regions of interest were computed. Notice the 
cluster of eigenvalues near the bottom edge of the plots: 
These are not resonances because they vary under pertur- 
bations in a. Instead, they belong to an approximation of 
the (scaled) continuous spectrum. 

It is more interesting to see log (N res ) as a function of 
— log (h) and R. This is shown in Figures [l9| and Us- 
ing least-squares regression, we can extract approximate 
slopes for the curves in Figure [l9|; these are shown in Ta- 
ble H and plotted in Figure [zi| 



■0.01 



0.02 



0.03 



0.04 



hbar=0. 022702 



O 



9j 



2i ^ 



0.3 0.4 0.5 0.6 0.7 



6.2 Trapped Set Dimension 

For comparison, £i£jp+l j s plotted as a function of R 
in Figure ^2| The figure contains curves corresponding 
to different energies E: The top curve corresponds to 
E = 0.4, the middle curve E — 0.5, and the bottom curve 
E = 0.6. It also contains curves corresponding to dif- 



Figure 15: Eigenvalues for E = 0.5, fi = 0.022702, R = 
1.4, and a G {0.0567, 0.0741, 0.0916}, using 112 x 119 
grid and 98 out of N = 13328 eigenvalues. 
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hbar=0.020616 
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Figure 16: Eigenvalues for E = 0.5, % = 0.020616, R = 
1.4, and a e {0.0515, 0.0689, 0.0864}, using 123 x 131 
grid and 107 out of iV = 16113 eigenvalues. 

hbar=0.018721 



® : * 





0.5 0.6 



Figure 18: Eigenvalues for E = 0.5, fi = 0.017, R = 1.4, 
and a G {0.0425,0.0599,0.0774}, using 149 x 159 grid 
and 127 out of N = 23691 eigenvalues. 



0.3 0.4 



Figure 17: Eigenvalues for E = 0.5, Ti = 0.018721, R = 
1.4, and a <E {0.0468,0.0642,0.0817}, using 135 x 144 
grid and 116 out of N = 19440. 
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E = 0.4 


E = 0.5 


E = 0.6 
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1.2558 
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1.3007 
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1.2521 
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1.2488 
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1.2636 


1.2524 




Figure 19: log(iV res ) as a function of — log(?i), for % 
varying from 0.017 to 0.025 and R = 1.4 + 0.05 • k, with 
< k < 6. (The lowest curve corresponds to R = 1.7, 
while the highest curve corresponds to R = 1.4.) 



Table 1 : Estimates of D ( K ^ +1 as a function of R. 
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Figure 22: This figure shows as a function of 

R: The top group of curves have E = 0.4, the middle 
E = 0.5, and the bottom E = 0.6. Solid curves marked 
with circles represent computations where No — 10000, 



Ni = 100, ^ < < % and 



< 



< 



Dashed 



curves marked with X's represent computations where 
A^o = 14142, whereas dashed curves marked with tri- 
angles represent computations where iVo = 10000 and 
Ni = 71. The recursion depth fco in all these figues is 6. 
The E = 0.6 curve does not appear to have completely 
converged but suffices for our purpose here. 



Figure 2 1 : The slopes extracted from Figure [T^, as a func- 
tion of R. The dotted curve is a least-squares regression 
of the "noisy" curve. 
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R 


slope 


D(K E ) + 1 
2 


relative error 


1.4 


1.2475 


1.2885 


0.032888 


1.45 


1.3433 


1.2834 


0.044645 


1.5 


1.2822 


1.2829 


0.00052244 


1.55 


1.327 


1.2773 


0.037472 


1.6 


1.3055 


1.2725 


0.025256 


1.65 


1.2304 


1.2694 


0.031756 


1.7 


1.2431 


1.2636 


0.016509 



Table 2: This table shows the slopes extracted from Figure 
|l9| , as well as the scaling exponents one would expect if 
the conjecture were true ({D (Ke) + 1) /2, computed at 
E = 0.5). Relative errors are also shown. 



ferent program parameters, to test the numerical conver- 
gence of dimension estimates. These curves were com- 
puted with 9 6 [^p, 4p], pg 6 [— | , |], and recursion 
depth ko — 6 (corresponding to symmetric sequences of 
length 2 ■ 6 + 1 = 13); the caption contains the values of 
Nq and Ni for each curve. For reference, Table [l] con- 
tains the dimension estimates shown in the graph. It is 
important to note that, while the dimension does depend 
on E and R, it only does so weakly: Relative to its value, 
d(Ke)+i j s ver y j-pygjjjy cons tant across the range of R 




Figure 23: Dashed lines with circles represent Eil£si+± 
as functions of R, for E e {0.4, 0.5, 0.6}. The dotted 
curve with triangles is the scaling exponent curve from 
Figure [2lJ while the solid curve with stars is the linear 
regression curve from that figure. Relative to the value of 
the dimension, the fluctuations are actually fairly small: 
See Table for a quantitative comparison. 



and E computed here. 



6.3 Discussion 



Table □ contains a comparison of 



D{K E ) + 1 



(for E = 0.5) 



R=1.4 



as a function of R, versus the scaling exponents from Fig- 
ure Figure 23 is a graphical representation of simi- 
lar information. This figure shows that even though the 
scaling curve in Figure |T] is noisy, its trend nevertheless 
agrees with the conjecture. Furthermore, the relative size 
of the fluctuations is small. At the present time, the source 
of the fluctuation is not known, but it is possibly due to the 
fact that the range of ft explored here is simply too large 
to exhibit the asymptotic behavior clearly |^| 

Figures ^4] - [30] contain plots of log (N res ) versus 
— log (Ti), for various values of R. Along with the nu- 
merical data, the least-squares linear fit and the scaling 
law predicted by the conjecture are also plotted]^] In con- 
trast with Figure ^ these show clear agrement between 
the asymptotic distribution of resonances and the scaling 
exponent predicted by the conjecture. 



17 But see Footnote 

18 The conjecture only supplies the exponents for power laws, not the 
constant factors. In the context of these logarithmic plots, this means the 
conjecture gives us only the slopes, not the vertical shifts. It was thus 
necessary to compute an ^-intercept for each "prediction" curve (for the 
scaling law predicted by the conjecture) using least squares. 




Figure 24: For R — 1.4: Triangles represent numerical 
data, circles least squares regression, and stars the slope 
predicted by the conjecture, h ranges from 0.025 down to 
0.017. 
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Figure 25: Same for R = 1.45. Again, T% ranges from 
0.025 to 0.017. 



Figure 27: R= 1.55 
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Figure 26: R = 1.5, 0.017 < h < 0.025. 



Figure 28: R = 1.6, 0.017 < h < 0.025. 
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Figure 3 1 : Resonances for two-bump scattering with Ti = 
0.035. 



Figure 29: R = 1.65, 0.017 < h < 0.025. 



R=1.7 
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Figure 30: R = 1.7, 0.017 < % < 0.025. 



6.4 Double Gaussian Scattering 

Finally, we compute resonances for the double gaussian 
model (setting m = 2 in (||). This case is interesting for 
two reasons: First, there exist rigorous results Jl3| , pT| | 
against which we can check the correctness of our results. 
Second, it helps determine the validity of semiclassical ar- 
guments for the values of % used in computing resonances 
for the triple gaussian model. 

The resonances are shown in Figures [H] - [37]: In these 
plots, R = 1.4 and h ranges from 0.035 to 0.015. One 
can observe apparent pseudospectral effects in the first 
few figures [ |30| , |3^ ]; this is most likely because the scal- 
ing angle a used here is twice as large as suggested in 
Section \. 1 , to exhibit the structure of resonances farther 



away from the real axis. 

To compare this information with known results [ pjj 
|IJ], we need some definitions: For a given energy < 
E < 1, define C(E) by 



C(E) = 2 1 v / 2-(E-V(x))dx, 

Jxo(E) 



where the limits of integration are 



x (E) 
Xi{E) 



= -R+ ^-2aHog(E ), 
= R-^/-2a 2 log(E). 



(48) 



(49) 



Let 0(E) denote the larger (in absolute value) eigenvalue 
of _D$(0, 0); log (6) is the Lyapunov exponent of $, and 
is easy to compute numerically in this case. Note that for 
two-bump scattering, each energy E determines a unique 
periodic trapped trajectory, and C (E) is the classical ac- 
tion computed along that trajectory. 
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hbar=0.030391 hbar=0.022913 




Figure 32: Resonances for two-bump scattering with h = Figure 34: Resonances for two-bump scattering with h = 
0.030391. 0.022913. 



hbar=0.026388 hbar=0.019895 




Figure 33: Resonances for two-bump scattering with h = Figure 35: Resonances for two-bump scattering with h = 
0.026388. 0.019895. 
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Figure 38: Lattice points for Ti = 0.035. 



Figure 36: Resonances for two-bump scattering with Ti 
0.017275. 



hbar=0.015 




Since these expressions are analytic, they have contin- 
uations to a neighborhood of the real line — C(E) be- 
comes a contour integral. In fll3|], it was shown that any 
resonance A = E — i-f must satisfy 



C(A) 



2ttH (m + ±) - 

ih (n+±)log (9 (Re (A))) 

o(n 2 ), 



(50) 



i in to + i 



where m and n are nonnegative integers. (The 
comes from the Maslov index associated with the classi- 
cal turning points.) This suggests that we define the map 

F (A) = (Fi (A) , F 2 (A)), where 



F x (A) 



and 



F 2 (A) 



ife (C(A)) _ 1 
2irh 2 

/to ((7(A)) 



(51) 



(52) 



7ilog(0(fle(A))) 2 
F should map resonances to points on the square integer 



contain 



0.3 0.4 0.5 0.6 0.7 



lattice, and this is indeed the case: Figures |38|- 44 ( 
images of resonances under F, with circles marking the 
nearest lattice points. The agreement is quite good, in 
view of the fact that we neglected terms of order h 2 in 
Equation (p0|). 



7 Conclusions 

Figure 37: Resonances for two-bump scattering with Ti = 

0.015. Using standard numerical techniqes, one can compute a 

sufficiently large number of resonances for the triple gaus- 
sian system to verify their asymptotic distribution in the 
semiclassical limit Ti — > 0. This, combined with effective 
estimates of the fractal dimension of the classical trapped 
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Figure 39: Lattice points for h = 0.030391. 
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Figure 41: Lattice points for h = 0.022913. 
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Figure 40: Lattice points for h = 0.026388. 
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Figure 42: Lattice points for ft, = 0.019895. 
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Figure 43: Lattice points for h = 0.017275. 



set, gives strong evidence that the number of resonances 
N res in a box [Eq,Ei] — i [0,h], for sufficiently small 
\Et - E Q \ andft, 



(53) 



as one can see in Figure ^3] and Table |5| Furthermore, the 
same techniques, when applied to double gaussian scatter- 
ing, produce results which agree with rigorous semiclassi- 
cal results. This supports the correctness of our algorithms 
and the validity of semiclassical arguments for the range 
of Ti explored in the triple gaussian model. The compu- 
tation also hints at more detailed structures in the distri- 
bution of resonances: In Figures 14 - [TJ| one can clearly 
see gaps and strips in the distribution of resonances. A 
complete understanding of this structure requires further 
investigation. 

While we do not have rigorous error bounds for the di- 
mension estimates, the numerical results are convincing. 
It seems, then, that the primary cause for our failure to ob- 
serve the conjecture in a "clean" way is partly due to the 
size of h: If one could study resonances at much smaller 
values of h, the asymptotics may become more clear. 
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